/*
 * Copyright (C) 2010,2011  Paulo Cesar Pereira de Andrade.
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * Authors:
 *	Paulo Cesar Pereira de Andrade
 */

#include "exl.h"

#include <math.h>
#include <signal.h>

/*
 * Types
 */
typedef struct {
    evector_t	*r;
    evector_t	*a;
    evector_t	*b;
    etype_t	 t;
} evop_t;

#if HAVE_SIMD
typedef eint8_t			vc_t __attribute__((vector_size(16)));
typedef eint16_t		vs_t __attribute__((vector_size(16)));
typedef eint32_t		vi_t __attribute__((vector_size(16)));
typedef eint64_t		vl_t __attribute__((vector_size(16)));
typedef efloat32_t		vf_t __attribute__((vector_size(16)));
typedef efloat64_t		vd_t __attribute__((vector_size(16)));
#  define vc_length(l)		((l + 0xf) >> 4)
#  define vs_length(l)		((l + 0x7) >> 3)
#  define vi_length(l)		((l + 0x3) >> 2)
#  define vl_length(l)		((l + 0x1) >> 1)
#else
#  define vc_t			eint8_t
#  define vs_t			eint16_t
#  define vi_t			eint32_t
#  define vl_t			eint64_t
#  define vf_t			efloat32_t
#  define vd_t			efloat64_t
#  define vc_length(l)		(l)
#  define vs_length		vc_length
#  define vi_length		vc_length
#  define vl_length		vc_length
#endif

#define vo_t			eobject_t
#define c_t			eint8_t
#define s_t			eint16_t
#define i_t			eint32_t
#define l_t			eint64_t
#define f_t			efloat32_t
#define d_t			efloat64_t
#define o_t			eobject_t

#define vf_length		vi_length
#define vd_length		vl_length

/*
 * Prototypes
 */
#define PROTO			1
#include "evector_set.c"
#include "evector_and.c"
#include "evector_or.c"
#include "evector_xor.c"
#include "evector_com.c"
#include "evector_add.c"
#include "evector_sub.c"
#include "evector_mul.c"
#include "evector_div.c"
#include "evector_trunc.c"
#include "evector_rem.c"
#include "evector_neg.c"
#include "evector_mul2.c"
#include "evector_div2.c"
#include "evector_shl.c"
#include "evector_shr.c"
#include "evector_signbit.c"
#include "evector_signum.c"
#include "evector_rational.c"
#include "evector_num.c"
#include "evector_den.c"
#include "evector_real.c"
#include "evector_imag.c"
#include "evector_conj.c"
#include "evector_arg.c"
#include "evector_abs.c"
#include "evector_sqrt.c"
#include "evector_cbrt.c"
#include "evector_sin.c"
#include "evector_cos.c"
#include "evector_tan.c"
#include "evector_asin.c"
#include "evector_acos.c"
#include "evector_atan.c"
#include "evector_sinh.c"
#include "evector_cosh.c"
#include "evector_tanh.c"
#include "evector_asinh.c"
#include "evector_acosh.c"
#include "evector_atanh.c"
#include "evector_proj.c"
#include "evector_exp.c"
#include "evector_log.c"
#include "evector_log2.c"
#include "evector_log10.c"
#include "evector_atan2.c"
#include "evector_pow.c"
#include "evector_hypot.c"
#undef PROTO

static void
setup(evop_t *o, evector_t *a, evector_t *b, etype_t at, etype_t bt,
      ebool_t i, ebool_t f);

static etype_t
match(etype_t a, etype_t b, ebool_t i, ebool_t f);

static evector_t *
create(evector_t *v, etype_t t);

static evector_t *
coerce(evector_t *v, etype_t f, etype_t t);

static void
set_pointer_int(eobject_t *u, eint_t i);

static void
set_pointer_float(eobject_t *u, efloat_t d);

static void
set_pointer_mpz(eobject_t *u, empz_t z);

static void
set_pointer_mpz_check(eobject_t *u, empz_t z);

static void
set_pointer_mpq(eobject_t *u, empq_t q);

static void
set_pointer_mpq_check(eobject_t *u, empq_t q);

static void
set_pointer_mpr(eobject_t *u, empr_t r);

static void
set_pointer_cdd(eobject_t *u, ecdd_t dd);

static void
set_pointer_cqq(eobject_t *u, ecqq_t qq);

static void
set_pointer_mpc(eobject_t *u, empc_t cc);

#define set_i(i)			set_pointer_int(u, i)
#define set_d(d)			set_pointer_float(u, d)
#define set_z(z)			set_pointer_mpz(u, z)
#define set_z_check(z)			set_pointer_mpz_check(u, z)
#define set_q(q)			set_pointer_mpq(u, q)
#define set_q_check(q)			set_pointer_mpq_check(u, q)
#define set_r(r)			set_pointer_mpr(u, r)
#define set_dd(dd)			set_pointer_cdd(u, dd)
#define set_qq(qq)			set_pointer_cqq(u, qq)
#define set_cc(cc)			set_pointer_mpc(u, cc)

#define get_i(x)			(**(eint_t **)x)
#define get_d(x)			(**(efloat_t **)x)
#define get_z(x)			(*(empz_t *)x)
#define get_q(x)			(*(empq_t *)x)
#define get_r(x)			(*(empr_t *)x)
#define get_dd(x)			(**(ecdd_t **)x)
#define get_qq(x)			(*(ecqq_t *)x)
#define get_cc(x)			(*(empc_t *)x)

/*
 * Implementation
 */
#define CODE			1
#include "evector_set.c"
#include "evector_and.c"
#include "evector_or.c"
#include "evector_xor.c"
#include "evector_com.c"
#include "evector_add.c"
#include "evector_sub.c"
#include "evector_mul.c"
#include "evector_div.c"
#include "evector_trunc.c"
#include "evector_rem.c"
#include "evector_neg.c"
#include "evector_mul2.c"
#include "evector_div2.c"
#include "evector_shl.c"
#include "evector_shr.c"
#include "evector_signbit.c"
#include "evector_signum.c"
#include "evector_rational.c"
#include "evector_num.c"
#include "evector_den.c"
#include "evector_real.c"
#include "evector_imag.c"
#include "evector_conj.c"
#include "evector_arg.c"
#include "evector_abs.c"
#include "evector_sqrt.c"
#include "evector_cbrt.c"
#include "evector_sin.c"
#include "evector_cos.c"
#include "evector_tan.c"
#include "evector_asin.c"
#include "evector_acos.c"
#include "evector_atan.c"
#include "evector_sinh.c"
#include "evector_cosh.c"
#include "evector_tanh.c"
#include "evector_asinh.c"
#include "evector_acosh.c"
#include "evector_atanh.c"
#include "evector_proj.c"
#include "evector_exp.c"
#include "evector_log.c"
#include "evector_log2.c"
#include "evector_log10.c"
#include "evector_atan2.c"
#include "evector_pow.c"
#include "evector_hypot.c"
#undef CODE

static void
setup(evop_t *o, evector_t *a, evector_t *b, etype_t at, etype_t bt,
      ebool_t i, ebool_t f)
{
    if (a->rank != b->rank ||
	memcmp(&a->dims[0], &b->dims[0], a->rank * sizeof(edim_t)))
	eraise_exception(except_invalid_argument);

    o->a = a;
    o->b = b;

    v_check(3);

    o->t = match(at, bt, i, f);
    if (at != o->t)
	o->a = coerce(o->a, at, o->t);
    if (bt != o->t)
	o->b = coerce(o->b, bt, o->t);
    o->r = create(o->a, o->t);
}

static etype_t
match(etype_t a, etype_t b, ebool_t i, ebool_t f)
{
    etype_t	t;

    /* if types don't match or an unsigned type, calculate as a wider
    * signed type if available, otherwise, resort to dynamic type;
    * typed vectors also means overflows are ignored, so, if types
    * matches, calculations are not done in a wider type */
    if (i) {
	if (a == t_float32 || a == t_float64 ||
	    b == t_float32 || b == t_float64)
	    eraise_exception(except_not_an_integer);
    }
    switch (a) {
	case t_int8:
	    switch (b) {
		case t_int8:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int8;
		    break;
		case t_uint8:	case t_int16:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int16;
		    break;
		case t_uint16:	case t_int32:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int32;
		    break;
		case t_uint32:	case t_int64:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int64;
		    break;
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_uint8:		case t_int16:
	    switch (b) {
		case t_int8:	case t_uint8:	case t_int16:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int16;
		    break;
		case t_uint16:	case t_int32:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int32;
		    break;
		case t_uint32:	case t_int64:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int64;
		    break;
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_uint16:		case t_int32:
	    switch (b) {
		case t_int8:	case t_uint8:	case t_int16:	case t_uint16:
		case t_int32:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int32;
		    break;
		case t_uint32:	case t_int64:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int64;
		    break;
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_uint32:		case t_int64:
	    switch (b) {
		case t_int8:	case t_uint8:	case t_int16:	case t_uint16:
		case t_int32:	case t_uint32:	case t_int64:
		    t = f ? !tmain->flt ? t_float64 : t_void : t_int64;
		    break;
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_uint64:
	    switch (b) {
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_float32:
	    switch (b) {
		case t_int8:	case t_uint8:	case t_int16:	case t_uint16:
		case t_int32:	case t_uint32:	case t_int64:	case t_uint64:
		case t_float32:
		    t = t_float32;
		    break;
		case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	case t_float64:
	    switch (b) {
		case t_int8:	case t_uint8:	case t_int16:	case t_uint16:
		case t_int32:	case t_uint32:	case t_int64:	case t_uint64:
		case t_float32:	case t_float64:
		    t = t_float64;
		    break;
		default:
		    t = t_void;
		    break;
	    }
	    break;
	default:
	    t = t_void;
	    break;
    }

    return (t);
}

static evector_t *
create(evector_t *v, etype_t t)
{
    evector_t	*r;

    enew_vector((evector_t **)&thread_self->obj, t, v->length, v->rank);
    r = thread_self->obj;
    v_push(r);

    if (v->rank != 1)
	memcpy(&r->dims, &v->dims, v->rank << 3);

    return (r);
}

static evector_t *
coerce(evector_t *v, etype_t f, etype_t t)
{
    eint32_t	 l;
    evector_t	*r;

    l = v->length;
    r = create(v, t);

    switch (t) {
	case t_int16:
	    switch (f) {
		case t_int8:	set_vi16_vi8(r->v.obj, v->v.obj, l);	break;
		default:	set_vi16_vu8(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
	case t_int32:
	    switch (f) {
		case t_int8:	set_vi32_vi8(r->v.obj, v->v.obj, l);	break;
		case t_uint8:	set_vi32_vu8(r->v.obj, v->v.obj, l);	break;
		case t_int16:	set_vi32_vi16(r->v.obj, v->v.obj, l);	break;
		default:	set_vi32_vu16(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
	case t_int64:
	    switch (f) {
		case t_int8:	set_vi64_vi8(r->v.obj, v->v.obj, l);	break;
		case t_uint8:	set_vi64_vu8(r->v.obj, v->v.obj, l);	break;
		case t_int16:	set_vi64_vi16(r->v.obj, v->v.obj, l);	break;
		case t_uint16:	set_vi64_vu16(r->v.obj, v->v.obj, l);	break;
		case t_int32:	set_vi64_vi32(r->v.obj, v->v.obj, l);	break;
		default:	set_vi64_vu32(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
	case t_float32:
	    switch (f) {
		case t_int8:	set_vf32_vi8(r->v.obj, v->v.obj, l);	break;
		case t_uint8:	set_vf32_vu8(r->v.obj, v->v.obj, l);	break;
		case t_int16:	set_vf32_vi16(r->v.obj, v->v.obj, l);	break;
		case t_uint16:	set_vf32_vu16(r->v.obj, v->v.obj, l);	break;
		case t_int32:	set_vf32_vi32(r->v.obj, v->v.obj, l);	break;
		case t_uint32:	set_vf32_vu32(r->v.obj, v->v.obj, l);	break;
		case t_int64:	set_vf32_vi64(r->v.obj, v->v.obj, l);	break;
		default:	set_vf32_vu64(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
	case t_float64:
	    switch (f) {
		case t_int8:	set_vf64_vi8(r->v.obj, v->v.obj, l);	break;
		case t_uint8:	set_vf64_vu8(r->v.obj, v->v.obj, l);	break;
		case t_int16:	set_vf64_vi16(r->v.obj, v->v.obj, l);	break;
		case t_uint16:	set_vf64_vu16(r->v.obj, v->v.obj, l);	break;
		case t_int32:	set_vf64_vi32(r->v.obj, v->v.obj, l);	break;
		case t_uint32:	set_vf64_vu32(r->v.obj, v->v.obj, l);	break;
		case t_int64:	set_vf64_vi64(r->v.obj, v->v.obj, l);	break;
		case t_uint64:	set_vf64_vu64(r->v.obj, v->v.obj, l);	break;
		default:	set_vf64_vf32(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
	default:
	    switch (f) {
		case t_int8:	set_vo_vi8(r->v.obj, v->v.obj, l);	break;
		case t_uint8:	set_vo_vu8(r->v.obj, v->v.obj, l);	break;
		case t_int16:	set_vo_vi16(r->v.obj, v->v.obj, l);	break;
		case t_uint16:	set_vo_vu16(r->v.obj, v->v.obj, l);	break;
		case t_int32:	set_vo_vi32(r->v.obj, v->v.obj, l);	break;
		case t_uint32:	set_vo_vu32(r->v.obj, v->v.obj, l);	break;
		case t_int64:	set_vo_vi64(r->v.obj, v->v.obj, l);	break;
		case t_uint64:	set_vo_vu64(r->v.obj, v->v.obj, l);	break;
		case t_float32:	set_vo_vf32(r->v.obj, v->v.obj, l);	break;
		case t_float64:	set_vo_vf64(r->v.obj, v->v.obj, l);	break;
		default:	set_vo_vo(r->v.obj, v->v.obj, l);	break;
	    }
	    break;
    }

    return (r);
}

static void
set_pointer_int(eobject_t *u, eint_t i)
{
    if (*u && eobject_type(*u) == t_int)
	**(eint_t **)u = i;
    else
	enew_integer(u, i);
}

static void
set_pointer_float(eobject_t *u, efloat_t d)
{
    if (*u && eobject_type(*u) == t_float)
	**(efloat_t **)u = d;
    else
	enew_float(u, d);
}

static void
set_pointer_mpz(eobject_t *u, empz_t z)
{
    if (*u && eobject_type(*u) == t_mpz)
	mpz_set(*(empz_t *)u, z);
    else
	enew_mpz(u, z);
}

static void
set_pointer_mpz_check(eobject_t *u, empz_t z)
{
    if (empz_fit_i(z))
	set_i(empz_get_i(z));
    else if (*u && eobject_type(*u) == t_mpz)
	mpz_set(*(empz_t *)u, z);
    else
	enew_mpz(u, z);
}

static void
set_pointer_mpq(eobject_t *u, empq_t q)
{
    if (*u && eobject_type(*u) == t_mpq)
	mpq_set(*(empq_t *)u, q);
    else
	enew_mpq(u, q);
}

static void
set_pointer_mpq_check(eobject_t *u, empq_t q)
{
    if (mpz_cmp_ui(mpq_denref(q), 1) == 0)
	set_pointer_mpz_check(u, mpq_numref(q));
    else if (*u && eobject_type(*u) == t_mpq)
	mpq_set(*(empq_t *)u, q);
    else
	enew_mpq(u, q);
}

static void
set_pointer_mpr(eobject_t *u, empr_t r)
{
    if (*u && eobject_type(*u) == t_mpr)
	mpfr_set(*(empr_t *)u, r, thr_rnd);
    else
	enew_mpr(u, r);
}

static void
set_pointer_cdd(eobject_t *u, ecdd_t dd)
{
    if (imag(dd) == 0.0)
	set_d(real(dd));
    else if (*u && eobject_type(*u) == t_cdd)
	**(ecdd_t **)u = dd;
    else
	enew_cdd(u, dd);
}

static void
set_pointer_cqq(eobject_t *u, ecqq_t qq)
{
    if (mpq_sgn(cqq_imagref(qq)) == 0) {
	if (mpz_cmp_ui(mpq_denref(cqq_realref(qq)), 1))
	    set_q(cqq_realref(qq));
	else
	    set_z_check(mpq_numref(cqq_realref(qq)));
    }
    else if (*u && eobject_type(*u) == t_cqq)
	cqq_set(*(ecqq_t *)u, qq);
    else
	enew_cqq(u, qq);
}

static void
set_pointer_mpc(eobject_t *u, empc_t cc)
{
    if (mpfr_zero_p(mpc_imagref(cc)))
	set_r(mpc_realref(cc));
    else if (*u && eobject_type(*u) == t_mpc)
	mpc_set(*(empc_t *)u, cc, thr_rndc);
    else
	enew_mpc(u, cc);
}
